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Abstract 

This introduction to Green's functions is based on their role as kernels of differential equations. 
The procedures to construct solutions to a differential equation with an external source or with an 
inhomogeneity term are put together to derive the Dyson equation for the Green's function of the 
inhomogeneous system. Very different areas of physics such as, for example, electrodynamics and 
quantum transport, can profit from this Green's function formalism. 

The fundamental homogeneous-medium Green's tensor of electrodynamics is deduced from the 
field of a dipole. Based upon that a numerical procedure is presented to solve the wave-equation 
for the near-field in a scattering setup for arbitrary material distributions. The full inhomogeneous 
system's Green's function is not explicitly needed to get the fields, although it can be obtained by 
a very similar calculation and in optics can be interpreted as a density of states. 

It is demonstrated how the transport problem for two open free-electron gas reservoirs with 
arbitrary coupling can be solved by finding the system's Green's function. In this sense the article 
is an introduction on Green's functions for treating interaction. A very detailed discussion of the 
current formula is given on an elementary basis. 
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I. GREEN'S FUNCTIONS TOOL FOR SOLVING DIFFERENTIAL EQUATIONS 



Green's functions [H [21 El S] are encountered as response functions, time-ordered ex- 
pectation values, certain solutions of boundary- value problems or resolvent kernels. This 
introduction to Green's functions is based on their role as kernels of differential equations. 
The procedures to construct solutions to a differential equation with an external source or 
with an inhomogeneity term are put together to derive the Dyson equation for the Green's 
function of the inhomogeneous system. Very different areas of physics such as, for example, 
electrodynamics (see section II and [3]) and quantum transport (see section III and [SIE]), 
can profit from such Green's function formalisms. 

1.1. Introduction 

Green's function formalisms do not present a cure-all for solving differential equations, 
because essentially the problem of finding the solution of the differential equation is shifted 
to that of finding the corresponding Green's function. This can, however, be a simplification 
and even give access to the solution of a more general class of problems. Starting from or- 
dinary differential equations, this short review is meant to introduce how the corresponding 
Green's functions are defined and how they are involved in constructing the solutions for 
different types of differential equations. Although bearing similarities, the Green's function 
formalism can go beyond perturbation theory. Special emphasize is laid here upon the par- 
allels between homogeneous and inhomogeneous ordinary differential and Green's functions 
equations. 

1.2. Homogeneous equation 

The starting point is a homogeneous differential equation 

0o(O = (1) 

which we suppose is exactly solvable, although 0o(O will not explicitly be needed. is 
some differential operator which may include multiplication with a constant or even another 
function of ^. ,^ is either a space or time variable. What we need to know is the solution g 
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of the corresponding Green's functions equation 

Dig{^,0 = S{^-0 (2) 

g is a tensor-like function of two arguments, only acting on the first of them. Like g 
replaces on the left, the zero on the right side of (1) is replaced by a ^-distribution in (2). 
There is no general recipe, but knowing 0o(O can help to get g{^,^'). 

1.3. Source term 

Having g, the construction of a solution 

+ j di'g{i,OQ{0 (3) 

of the differential equation with a source term Q on the right side 

MO = QiO (4) 

is straight forward. Of course, a solution 0o of (1) can be added independently of Q, so we 
only need to proove that the integral term from (3) satisfies (4): 

= J d^'6{^-0 Q{0 = QiO 

(All integrals are understood to range over the entire ^-space.) 

1.4. Inhomogeneity 

Instead of a source term, the differential equation can contain a potential term which we 
shall call an inhomogeneity. 

MO - viO MO = (5) 

that is a ^-dependent function V multiplied with 0. V{^) is not included in D^, because 
we assume that it so much complicates the equation that a standard solution is no longer 
known. The minus sign is a useful convention. We put the inhomogeneity term on the right 

MO = V{0 MO (5a) 
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to make (5) formally look like (4). With V(f) playing the role of Q a formal solution is 
constructed analogously to (3): 

MO = MO + J dO g{i, 0) v{0) MO) (6) 

This presents an implicit equation for (pi^, the so-called Lippmann-Schwinger equation. We 
have replaced the differential equation by an integral equation. Inserting the solution of 
the homogeneous equation 0o also for (pih in the integral on the right side of (6) would give 
the Born approximation 4>ih{0 ~ 0o(O + / ^0 diOO) ^{0) 4>o{0)^ which is appropriate 
if V(j) is a small perturbation compared to D^(f). However, it is the virtue of the Green's 
functions method that in contrast to perturbation theory the inhomogeneity need not be a 
small deviation. Although not yet providing an explicit solution for (pih, one can make use 
of (6) in numerical calculations (see section II and [8J). As (2) is the corresponding Green's 
function equation with a point source to the homogeneous equation (1), the equation defining 
the Green's function G for the inhomogeneous case derived from (5) reads: 

[Di-V{O]G{^,0) = S{^-0) (7) 

However, while supposing that we know or can easily guess g{0 0)y there is no hint yet how 
to calculate G{^,0)- 

1.5. General case 

Finally we have to treat the most general case with source and inhomogeneity term: 

MiO - viO MO = QiO (8) 

In analogy to (3) and (4) a yet formal solution can be written down using (7): 

MO = MO + J dO G{^, 0) QiO) (9) 

Here we have the freedom to add any solution of (6) to the integral. 

But let us start again solving the problem directly from (8) which we rewrite as 

MiO - Q{0 = viO MiO (8a) 
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We compare this with (5a) and construct a solution of the same form as (6). 0g takes the 
role of because it would satisfy (8a) if the right side were zero. has to be replaced 
by (t>ihq- Using (3) for (pq we obtain 



(10) could also have been established by a different hne of thought. In (8) putting everything 
that differs from the homogeneous equation (1) on the right side gives 

If it were not for the inhomegeneity term V0, the solution would be (3). And if it were 
not for the source term Q, we could use (6). Although this is, of course, not a correct 
way to solve non- homogeneous differential equations, we can understand (10) as an ansatz 
adding these two contributions. The ^o-part does not have to be written twice. And (pihq 
has to appear instead of (pit in the integral taken from (6). This still unknown (pihq leaves 
the neccessary freedom to somehow counterbalance the Q-term not present in (5) and (6), 
which justifies (10) as an ansatz. Whichever way it was obtained, (10) is an implicit integral 
equation for (pi^q as (6) is for (p^h- 

1.6. Dyson equation 

Now we shall profit from the fact that with (9) we have a second representation of (pihg. 
Insert the expression from (9) for (pi^q both on the left and on the right side in (10): 

MO + J di' G(e,e') Q{0 = MO + J dO 9(^,0 QiO + 

I degiW) vio') [ MO + / do G(r,e') Qio i 

(pih on the left cancels with 0o and one integral term on the right according to (6) and we 
are left with 

dO G{U') Q{0 = / dO g{i,0) Q{0 + 

di" j dOgii,e) viC) Gie,e)Qie) 
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This is an implicit equation for the unknown Green's function G. But as G by definition (7) 
does not depend on any source term, (11) must be vahd for arbitrary Q. The special choice 



(11) is the general implicit integral equation for the inhomogeneous system's Green's function 
G and called Dyson equation. 

We could have derived (11) by doing all calculations from (8) on with only a point source 
6{C, — ^o) instead of Q{^). Putting in the solutions (3) and (9) for this case would just have 
looked a bit awkward. (11) is not an explicit solution for G, but still an implicit equation, 
and even though it is often written as G = g + gVG one must not forget that there is a 
convolution-like integration over the inner ^-index behind the sequence of factors gVG. The 
integration can either be transformed into a discrete finite sum, and thus (11) into a linear 
system of equations solvable by a matrix inversion (see section II and [9j), or the convolution 
can be replaced by a multiplication going to frequency space by a Fourier transformation 
(see section III and 0). 

(10) and the preparation for (11) would have looked much more elegant leaving out (po 
and the 0j/i-terms right from the start. Indeed, one could argue that in (3) one is only 
interested in the part different from the trivial homogenous solution, if no further boundary 
conditions have to be accounted for. With a similar argument, the 0i/i-part could have been 
dropped in (9). The identities 



= - Co) gives 




which after carrying out the 



integrals and then renaming to ^' becomes 




(11) 




(3*) 




(6*) 




(9*) 



6 



(10*) 



would obviously also have given (11). Leaving out 0o in (6) would be precarious (see 
below). However, dropping it when constructing (11) along our second line of thought 
causes no problem. The term with Q will ensure that 0*^^ cannot simply be zero. 

1.7. Expansions 

There is some subtlety about the 0o contribution in (6). (6*) has the trivial solution 
4'ihiO — not necessarily another one. Surely, (f>ih{0 = satisfies (5), but it is not 
what we want. 0o is the background excitation replacing boundary conditions in this kind of 
problem, which shall become clear when discussing electrodynamics (see following article). 
00 acts as a source term which - instead of being put in Q - can be more appropriately and 
conveniently set as a fixed part of the we are looking for. Inserting (6) into itself ever and 
ever again, (pih is developed into a series in powers of gV: 



Just for shorthand notation we dropped the arguments and integration variables and in the 
third line also the integral signs. The resulting series is recognized from perturbation theory 
summing interactions to zeroth, first, second, etc. order. If there were no contribution (pg in 
(f>ih, only the contribution of gV to power infinity would exist with nothing to multiply to 
at the end. In other words we would have no basis on which to develop (f)ih. The Green's 
functions formalism as we use it aims at solving (6), (10) or (11) in a closed form, not cutting 
the series in (12) at some finite order. The development (12) is only shown here precisely to 
demonstrate that interactions are included to all orders as well as to explain the importance 
of the homogeneous background contribution. The Dyson equation (11) can be expanded 
in an analogous manner to the Lippmann-Schwinger equation, by the way prooving the 




))) 



oo 




(12) 



n=0 
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equivalence to its complementary form G = g + GVg: 



oo oo 



G = g + 



gVG = gJ2iVgr = J2^9V) 



n 



9 = 9 + 



GVg 



(13) 



n=0 n=0 



(Integration over inner arguments is understood in all contributions to the sums.) Either 
formally or as a matrix calculation in finite discrete ^-space, (11) is often solved as 



Developing the (. . .)~^-factor in (14) into a geometric series just results in the infinite sums 
written in (13). 

Putting together (12) and (13) we can get an alternative representation to (6) for the 
solution of the inhomogeneous differential equation: 



Therefore in the case that should it be easier to get the Green's function G than to solve 
the implicit equation (6) for (pih, we see that G can be useful also for treating the equation 
without an external source Q. Nevertheless, (15) again illustrates that Vcpo plays the role 
of the source. 

1.8. Conclusions and Outlook for section I 

Using Green's functions it has been shown how differential equations can be treated 
that differ from easily solvable ones by an additional potential term or an arbitrary source 
term. The formalism as presented here is for open-boundary in contrast to boundary-value 
problems 0, ID]- From constructing the solutions of the differential equations, we also 
obtained the constituting relation for the inhomogeneous system's Green's function, which 
characterizes the response to a point source and includes all-order interactions. Examples 
of applications will be given in the two following sections. 




(14) 




(15) 
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II. NON-BOUNDARY VALUE PROBLEMS IN NEAR-FIELD OPTICS 



The fundamental homogeneous-medium Green's tensor of electrodynamics is deduced 
from the field of a dipole. Based upon that a numerical procedure is presented to solve the 
wave-equation for the near-field in a scattering setup for arbitrary material distributions. 
The full inhomogeneous system's Green's function is not explicitly needed to get the fields, 
although it can be obtained by a very similar calculation and in optics can be interpreted 
density of states. 

11. 1. Introduction 

The typical problem in nano-optics p!T| [12] is the situation that some tiny structures are 
illuminated by an extended source, a plane wave for example, and then one is interested in the 
field distribution that arises from multiple scattering [13], especially to identify places where 
the field intensity gets considerably enhanced [T5] . A theory can be based on Green's 
functions, however, their implication differs slightly from the standardly taught cases of fixed 
boundary field values |16] or located sources. Modern optical scanning microscopes make 
it possible to probe and map directly even different quantities of the near-field [17], such 
as the electric and magnetic field intensities [IH] or the density of states [IH]. Applications 
of tayloring nano-structures with respect to optical properties include resonant particles 
jSniET] and cavities [221 [23], squeezed fields [21], wave guides and their adressing [25 l [261 [27] 
as well as transmission apertures [231 EH] and lithography masks [26] . We here present the 
Green's functions formalism that forms the bases of a finite-element quite effective numerical 
algorithm used in current research [2S1 ISHl EI]- This treatise is also given as an application 
example of the general methods to solve differential equations with certain perturbations 
presented in the preceeding paper. 

11. 2. Problem 

The discussion can be reduced to monochromatic light, that is a single frequency u and 
thus time dependence e"""^* for the fields. With non-magnetic materials the problem is to 
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— * 

find the solution E(r) of the wave equation 




(1) is a differential equation of the type (5) from section I. There is no source term on the 
right hand side of (1). If the source were, for example, a dipole located at some point, a 
source term with its oscillation strength and direction would have to be put on the right as 
Q6{r — tq). However, we shall see that a plane-wave source can be and is better included 
in (1) as it is. For simplicity we shall assume that the background medium, in which 
objects with different permittivities e are located (Fig.l), is vacuum with permittivity e^. 
For another embedding medium, its dielectric constant Sb would take the role of Eq. In 
(1) e{r) of the material distribution designates the dimensionless relative permittivity with 
respect to vacuum or the background medium. To separate (1) into a part representing a 
homogeneous differential equation with known solution and an inhomogeneity write its as 

2 2 

- V X V X E{r) + % E{r) + % {e{r) - 1) E{r) = (la) 
The correspondances to the quantities of the general formalism given in section I are 

2 2 

e^f, D^...^-VxVx... + ^- ..., V ^ -"^{e - 1), (Pih ^ E 

and 00 will become Ef,. We can write down the solution following section I after having 
prepared the background Green's function in the next section. 

II. 3. Background Green's function 

In any case we need the Green's function g{r, r') of the homogeneous problem satisfying 

- V,^ X Vr? X g{r, ?) + k'^g{f, ?) = 1 S{f- ?) (2) 

with k = uj/c. g will be a 3 ® 3 tensor or matrix here. More commonly [32] the small letter 
g is used for the scalar function 

_^ ^ik\f—r'\ 

9scai{r,r') = -—— — — fulfilling 
47r|r — r'l 

A,^ gscai {r, r') + k'^gscai {r, r') = 6{r- r') (3) 
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FIG. 1: Objects distributed in space and discretized into cubes of equal size. A plane in space 
discretized into a quadratic mesh is also shown. Arrows indicate only some possible scattering 
paths to a location marked by the cross. 

with k = uj/c. The Green's function g from (2) is then named with index h for ho- 
mogeneous. There are several ways to obtain g. One is based on the knowledge that 
if we have a scalar function \E'(r) solving + A;^ = 0, then Fi = V x (a\E') and 
F2 = V X V X (a\E') with a constant but arbitrary pivot vector a will both solve the vectorial 
equation —V x V x F + k^F = 0. The tensor g looked for in (2) can be constructed out 
of Fi{r), Fi{r'), F2{r), F2{r') together with F^lr) = V\l'(r) and Fs{r'). We shall not enter 
into the details of this mathematically slightly precarious approach [33j. A second recipe 
just mentioned here for completeness is given by the following statement [32]: If gscaii'f^,^') 
satisfies (3), then 

g{f, ?) = (1 + ^^^^) gscai{r, P) (4) 

is the tensor defined by (2). Of course, because of being in homogeneous space gscai and 
g effectively are functions oi R = f — r' alone. 1 is the unit matrix in 3 by 3 cartesian 
coordinate space and 
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( d d d 



V ® V means building a matrix out of derivatives 
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We shall deduce g from a physical reasoning. From standard electrodynamics [3^ one 
has the electric field of an oscillating dipole 

^2gifcR f -^2 _ ^(^^ 3R{Rp} - pR^ 



E{R) 



AireoR 



i?2 



R^ 




(5) 



It is important to take the exact formula here including retardation in contrast to common 
near- or far-field approximations. (5) gives the space part, the time dependence is just e~*'^* 
everywhere. To get the Green's tensor, we have to evaluate from (5) what field components in 
X-, y- and z-direction a dipole p at r' oriented along x would produce at r, what components a 
dipole oriented along y would produce and what components a dipole along z would produce 
and assemble all these in a matrix. A point dipole is the elementary excitation corresponding 
to the S on the right side of (2). The physical meaning of g is to tell us what field any such 
dipole would have. That is shown formally in the first matrix in (5a). Decomposing any p 
into its cartesian components, (5) can be rewritten as 



E{R) 



Ey{r) 
EJf) 



( 



2AkR 



AtteoR 



1 - 



1-ikR 



Px[r'j 
Px{r') 
Px{r') 

^^1 0^ 

1 
yOOly 



Ey{r) 
EJr) 



Py 



Py{r') 
Pyir') 
Py{r') 



EJf) 



Eyif) 
EJr) 



Pz[r 
Pz{r') 
Pz{r') J 

( 



'^^ (p\ 



Py 
\Vz) 



-?> + ?>ikR + k'^R^ 



X2 XY XZ 
YX YZ 
ZX ZY Z^ 



Py 
\PzJ 



(5a) 



from which we easily see that the matrix to be multiplied with p to produce E{R) is 



1 



ikR\ -3 + 3ikR + k'^R^ 

- R(S)R — 



ATreoR \ \^ k'^R^ J k^R'^ 

R without arrow means the absolut value |-R| = \^X^ + Y'^ + Z'^ and X,Y, Z stand for x—x', 
y — y' and z — z\ respectively. R^ R is the matrix written with X, Y and Z from (5a). 
The terms from (5) vectorially oriented along p cause the diagonal matrix contribution, 
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those stemming from terms with R^Rp) the full matrix in (5a). Compared to the above 
expression g{R) from (4) has a minus sign and misses a factor k'^ /eq. As will be discussed 
later, the source to put into equation (1) corresponding to an oscillating dipole is not the 
dipole moment p itself, but — /iocj^ times p. And because ~f^~2 = — 1 we have 

g.fc/j / / i_^kR\ ^ ^ -3 + 3tkR + k^R^\ 
'^''^ = -4^V V-^^)-''^'' k^^ ) 

The formula (6) fails for f = r' oi R = 0. g{R) including the case R = can be represented 
using the principal volume method [32j . In practice, working with finite elements, the value 
to put for g of its two arguments the same place can be derived from the polarization of a 
dielectric body. The discussion of g{r, r) is postponed to the next section. 



II.4. Solution for the field 



The starting point to find a solution for the electric field with the objects present is eq. 
(6) of section I, which rewritten in the variables of our problem here reads 

E{r) = Eb{r) + j rfV g{f,?) V{?) E{?) (7) 

where Ebif) is a solution of — V x V x ^^(r) + ^^^(r) =0 or already assumed to be the 
space part of a linearly polarized plane wave -Efe(r) = Eqc^''^ with a fixed amplitude vector 
Eq. The time dependence e"**^* can be omitted in Eb{r) as well as in E{f). 

From (7) a numerical procedure can be deduced if the objects with e{f) ^ Eq only 
occupy fractions of space rather small on the scale of the wavelength, not at all principally 
necessarily much smaller than A, though. Then we divide them up into finite elements 
(Fig.l) of volume Av, which we enumerate i = 1,. . . ,N and associate permittivities s{fi) 
or perturbances V{ri) = —^{e{ri) — 1) and local fields E{ri) uniform over Av. The linear 
sizes of the object cells should not exceed about A/10. They need not at all be placed on 
a regular grid. And the only reason for demanding small enough objects is not to get too 
many elements N. The integrand in (7) only exists at places where V does not vanish, and 
to evaluate the field E{r) also at such places r G {rj}, no values outside the objects appear 
in the equation. Changing to finite elements we thus get a linear system of equations for 
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the fields in the object cells 

N 

E[n) = Eb{ri) + 5^ At; g{fi, f,) V{f,) E{fj) (8) 
i=i 

which can be solved by a matrix inversion. To evaluate the resulting field at any other place, 
that is outside the objects, the E{fi) just have to be inserted into the finite-element version 
of (7): 

N 

E{f) = Ef,{f) + 5^ At; g{f, fj) V{fj) E{f,) r ^ {fj (9) 

Of course, we already needed g{fi,rj) with ri = fj to set up the system (8). Let us 
suppose that there is only a single cell with an e differing from the background eo- This 
is placed into a homogeneous field E^. If the cell has the shape of a sphere, the local field 
throughout its inside is aligned in the direction of Ef, and its value is Eioc = 2+1-^^ E^l ES]- 
No retardation effects have to be considered here as the size of the cell can in principle 
be made arbitrarily small. In (8) only keeping the term of the sum with i = j, setting 
E{fi} = Eioc and V{ri) = -P{e - 1) gives 

^ E, = E,- Av g{n,n) k^e - 1) ^ E^ 

from which follows that 

^^''"''•^ = 3ik, ' 

The factor 1/3 is also valid for cubic elementary cells, however, other shapes require different 
depolarization factors [321 ES] ■ 

The background and the resulting field at each point are already 3-vectors. Nevertheless, 
in order to solve (8), imagine the Es for the object cells assembled into long or "double" 
vectors of times 3 components 

E with (E)i = E(fi) and Eb with (Eb)i = Eb(fi)- 

With further the big 3A^ x 3N matrix M consisting of 3x3 blocks 

Mij = 1 Sij - Av g{ri,rj) V{rj) 

(8) then reads 

ME = Eb or E = M^^Eb (8a) 
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M can be inverted using the procedure described in the appendix, however, as M is a 
full matrix and the complete invers is needed, that is of no advantage and a standard 
inversion algorithm will do as well. 

For the problems of a few small scatterers here by introducing finite elements the implicit 
integral equation (7) for the electric field has been turned into a linear system of equations 
that is easily solvable. Even shying this effort, the coarsest, so-called Born approximation 
consists in replacing E in the integral in (7) or in the sums in (8) or (9) by E},. Keeping the 
summation over finite elements to estimate the integral we directly get 

N 

E{r) ^ E,{r) + din, rj) V{fj) E,{rj) (H) 

No distinction between places r in or outside the object cells is necessary in (11). The Born 
approximation only takes into account first-order scattering off every object and thus can 
only be good for weak scatterers with distances between them rather large on the scale of the 
wavelength. Producing a clearly different field pattern from the exact solution including all 
scattering orders, Fig. 2 demonstrates that the Born approximation is likely to be insufficient 
to model near-field optics setups. 

II.5. System Green's function and density of states 

Adding an arbitrary source term to our original wave equation (1) changes it into 

xE(r) + —e(r)E(f)^Q(r). (12) 

If now we know a tensor function G{f, r') satisfying 

-Vr'X Vr- X G(f, r') + —£{r)G{r, r') = 1 5{f - r') or (13) 
-V^- X V^. X G{f, f') + k'^Gif, P) - V{r)G{f, P) = 16{f- P) 

then obviously 

E{r) = J d^P G{f,P) Q{P) (14) 

would give a special solution of (12). Any solution of (1) could be added. As generally 
deduced in section I the implicit relation to get G from is the Dyson equation 

G{r, P) = g{r, P) + J d^P' g{r, P') V{P') G{P', P) (15) 

15 
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FIG. 2: Three cubes of 25nm side length and e = — 16.0 + l.Oi in a background of e = 1 illuminated 
by a p-polarized plane wave of unit amplitude with A = 500nm from the front and under an angle 
of 30° from beneath the plane defined by the cubes. Setup (a) and electric field intensity \E\'^ 
calculated via the Born approximation (b) and the exact Green's functions method (c). Besides 
subtle differences in the near-field pattern remark that in (b) the grey-scale is from 0.978 to 1.022 
whereas in (c) it is from 0.995 to 1.005. The intensity map is taken at a height 25nm above the 
cubes and given for an area of 3/im x S/zm. 
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G{f, r') for both arguments r'and r' covering all space is too much information to display 
at once and usually much more than what one is interested in. The imaginary part of G{r, r) 
is proportional to the density of states p [311 ES] . The deduction of this statement found in 
quantum mechanics book [37], however, rather argues with a system of energy eigenstates 
and the variation of the Green's function as well as the density of states with energy. No real 
f^-space is explicitly mentioned. Our interest lies in the spatial dependence of the density of 
states at fixed light frequency u. 

Even if described in terms of fields, concepts like reactance and work known from electrical 
circuits may be applied [381 EH]. The time average of the work done by the fields is given by 

Re ^ y d^r' J*{?) E{?) (16) 

With no other imposed fields, charges or currents than an oscillating point dipole, the 
latter will present the only external current J, which will thus be located as 5{r' — r). 
If the dipole moment oscillates as p{r',t) = pQe~'^'^^6{r' — r), the corresponding current is 
J{r',t) = ~iupoe~'^'^^6{r' — f). Deducing the wave equation for time harmonic fields (in 
vacuum for simplicity here) from Maxwells equations with current term 

V X ^ — iujfioH = and V x hqH + sopoiujE = p^J 

leads to 

^ ^ ^ cu^ - 
- V X V X E + —E = -iujuoJ (17) 

from which we see that the source term for the dipole has to be set as Q{r') = —iujfioJ{r') = 
—uj'^fiQPQe~^^*6{r' — r). The integral (16) reduces to the value of J*E at r. The electric field 
we get from (14): 

E{r) = -c<;Voe"'"*y" rfV G{f,r') S{r' - r) Po = -c<;Voe~*'"*G(r, f)po 

Inserting J* and E into (16) the time factors cancel as expected for a time average and but 
for a factor u^po/2 we get 

Re i Po G{r, r) Po = — Im pq G{f, r) pq (18) 

Choosing unit vectors along the coordinate axis for the probe dipole po, (20) will filter 
out the trace elements of the matrix G{r,r). We associate Px{r) oc ~lmGxx{f^,r), Pyif) oc 
-lmGyy{f,r), pz{r) oc -ImG^^ (f,f) and a total p(f) = px{r) + Py{r) + p^(f). 
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A motivation for taking the negative imaginary part of G(r, as a measure for the 
presence of modes can also be obtained by comparison to the energy resonance of a forced 
oscillator [101 • For optimal excitation from the energy point of view - in contrast to amplitude 
resonance - the force has to be 7r/2 ahead of the elongation or in phase with the velocity 
of the oscillator. G{f,f) describes the field caused by backaction of the system at the place 
of the probe dipole moment (taken as reference phase zero), and therefore —lmG{'r,f) is 
the part that can in a resonant manner further enhance the dipole oscillation. (In reality 
radiation out of the system will provide strong damping.) 

We now intend to evaluate a map of G{f, r) on, for example, a horizontal plane. The plane 
may lie above or below object cells or even cut some. Like the objects, the, of course, finite 
area of interest on the plane is divided into cells. Just depending on the desired resolution 
of the map the unit cell length of this mesh may well differ from the cell size chosen to 
discretize the objects (Fig.l). The list of object cell midpoints {f.i} from the last section, 
which shall be called region A, is extended by all cell midpoints from the map in the plane, 
which shall be called region B and is now understood to to be included in counting i from 
1 to a new N . Analogously to (8) the integral in (15) is replaced by a sum: 

N 

G{fi, fk) = g{ru n) + ^Av g{ri, fj) V{fj) G{fj,rk) (19) 

i=i 

It does not matter that the map B has a different mesh from Av as for the objects A, as 
V{rj) = for fj in B, anyway. (Spatial overlap of cells from A and B and even coincidence 
of midpoints is no problem; a place can be counted with V{fj) in A and without in B.) 

Although only G{fi, fi) with fi in region B is wanted as a final result, (21) has to be 
set up as an equation for a matrix of all G with each of its arguments any cell in A or S, 



schematically sketched as 



AA 


AB 


BA 


BB 



. To solve (21) for G we have to invert the same 



kind of matrix as in (8), the only difference being that j now also runs over the plane cells 
in addition to the object cells. 

G = M-^ g (21a) 

G and g are themselves matrices on contrast to vectors E and Eh- G consists of 3x3-blocks 
G{fj,fk)., g is made of blocks g{fi,fk) and the 3x3-block at position (^,j) in M given by 
16ij — Av g{fi,fj) V{fj). The g in (21a) and M like the big G-matrix have the structure 
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AA 


AB 


BA 


BB 



One could invert M as given, for example by the procedure from appendix 
A. The matrix M = 1 — M there is initialized with gV. Its AB and BB quadrants are zero 

[aa o\ 

and will stay zero throughout the procedure, M = . This is no contradiction, 

\^BA J 

as it is not M that is singular. Quadrant BA will be needed for multiplication with g in 
(21a). However there is an even more efficient algorithm to get G that already includes the 
multiplication by g. It directly calculates G = (1 — gV)~^g, which is the compact way to 
write (21a) as the solution of (15), also denoted G = g + gVG for short. The technical 
details can be found in appendix B. 

Trace components of G(r*, r) meaning densities of states for the three polarization di- 
rections (Fig. 3) above an optical coral |3T] in analogy to a quantum coral [H] have been 
measured [19] in a so called forbidden-light near- field optical microscope [12]. The sample 
consists of a stadium arrangement of gold particles on a glass surface. The forbidden-light 
setup prevents detecting light emitted from the fiber tip that has not passed through sur- 
face modes that make up the density of states for this system. Like for antinodal and nodal 
points in a resonator, more energy can go into the system when the excitation is placed at 
a point of high density of states than when coupling is bad where the density is low. 



II. 6. Remarks on the source terms and alternative solutions 

In section 4 we saw that it is convenient to start from a solution for the field in the form 
(7) if the excitation comes, for example, from a background field belonging to a plane wave. 
Though the matrix to invert bore a certain similarity to the evaluation of the Green's tensor 
in section 5, with (8) and (9) we directly calculated the field. In contrast, more adapted to 
localized sources, there is (14) as a solution of (12). If there is no additional background 
field to cause any excitation, no solution of the equation (1) with zero right side is to be 
added as further contribution and (14) is the field distribution to be observed. (14) has to 
be rewritten in terms of finite elements in order to be used in a numerical calculation. In 
the same way as the objects the source Q has to be devided into discrete cells or elementary 
dipoles. To distinguish their L locations from those of the objects we shall enumerate them 
as Pa, a = 1, 2, . . . , L. The place r to evaluate the field E may be anywhere outside or inside 
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FIG. 3: (a) Top view of the structure: gold pads (90x90x30nm) on a glass subtrate. Calculated 
densities of states pxx (b) Pyy (c) and pzz (d) in a plane lOOnm above the substrate for wavelength 
A = 543nm [reproduced after [I9J]. The grey scale is from -50 to 35 in (b), -45 to 30 in (c), -7.5 to 
15 in (d) and higher or lower values are white or black, respectively. 

the objects as well as beside or even at a source location. For the following development the 
Dyson equation for the Green's tensor is needed in a discretized form for both its variants 
G = g + gVG and G = ^ + GVg. 

L 

E{r) = Y,G{r,Pc.)Q{Pa) (22a) 

a=l 

L L N 

= g{f, p.)Q{pa) + E E ^jMrMrj, Pa)Q{po) (22b) 

L L N 

= E Si^' MQipo) + E E (^i' MQ{po) + 

a=l a=l j=l 

L N N 

E E E r^)V{r^)G{f,, f^)V{fj)g{f^, p^)Q{pa) (22c) 

a=l j=l i=l 
L L N 

= E 3^^' MQipo) + E E 3^^' ^j)Vi^j)G{rj, pa)Q{pa) 
a=l a=l j=l 

L N 

= E 9{r, p^)Q{p^) + E 9ir, r,)V{f,)E{r,) (22d) 

a=l j=l 

Having in mind a region where and a resolution with which E{r) is to be evaluated like 
the discretized plane B from the last section, it would be possible to supply G for all needed 
combinations of arguments (r, pa) and calculate E{r) as the single sum from (22a). To weave 
in the influence of the objects, G would have to be set up as a big matrix like in the last 
section over all combinations of three regions A, S and B here, the objects, the source and 
the map. Having calculated G in the A-iJ-scheme from the last section, one could evaluate 
E{f) with G{r,fj) as written in (22b). However, the most efficient way is given in (22c). 
G{fi,fj) is merely needed for and fj from the set of object cells, keeping a matrix to be 
inverted as small as possible, namely of {AA)-tjpe. Choosing and in (21) in the object 
set A, instead of in A or i? as the equation was originally set up for, we see that (21) presents 
a closed system of equations for all such G{ri,rk)- Then for (22c) more summations over 
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products with g^-functions, which are analytically known for any pair of arguments, can be 
considered less demanding in computing time than the inversion of large matrices. 

In the transformation from (22c) to (22d) after swapping index names i and j in the last 
sum, G = g + GVg and (14) have been exploited. Using (22d) for E{r) renders an implicit 
equation for the field in the form (10) or (10*) from section 1. E{r) is the equivalent of (pihq 
and as stated earlier, we assume that physically there is no background field 0o that could 
initiate an additional field distribution (p^h- Although not very convenient, a plane wave as 
exciting field could be understood as stemming from a sufficiently long and dense array of 

— * _ 

Huygens elementary dipole sources Q reasonably far away from the objects. The other way 

round, for a single dipole source or a number of dipole sources distributed in space the field 
they would produce at any location r in homogeneous space is the superposition of their 
individual fields, and putting Eh{f) — '^a=i di^^ PodQiM ansatz (7) can be used also 
for this case. 

— * 

Like the field E anywhere was obtained as a straight-forward summation once having its 
values at the places of the object cells, finally an alternative way to the procedure from the 
last section to get the Green's tensor G shall be given, also requiring only the inversion of a 
matrix with size the number of object cells. Series expansion is used to rewrite the solution 
of the Dyson equation: 

oo oo 

G = (1 - gVy'g = [1 + J2(9Vr]9 = 9 + ^^(^ W")^ 

n=l n=0 

= g + gVil-gV)-^g (23a) 

oo 

- 9 + 9Vg + gV{Y^{gVr)gVg = g + gVg + gV{l- gVy'g Vg 

n=0 

= g + gVg + gVGVg (23b) 
Designating regions the spatial arguments belong to on (23a) we get for Gbb- 

Gbb = gBB + gBA Va (1 - gV)']^A 9ab, (23a') 

As V does not vanish only in region A, the first index of (1 — gV)~^ obviously must be A. 
(1 — gV)~^ being X^^o(5'^)"' power zero the second region index automatically is the 
same as the first and all other powers ending with V imply second index A. It is sufficient 
to set up the matrix 1 — gV as an (7474)-block and invert that. Should one prefer to evaluate 
a complete G, which differs from the above inverted matrix by a factor g, line (23b) like 
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(22c) shows that it is in principle only necessary to get G from some self-consistent implicit 
equation in the object region A. Gaa can be constructed applying the procedure described 
in appendix B to a matrix set up as (AA)-block only. The inversion has to be completed 
in this case, though. Going through the diagonal elements, all lines and columns have to 
be updated in each step, including the ones above and to the left of as well as the ones the 
respective diagonal element is in. Writing (23b) 

Gbb = gsB + Qba Va Qab + Qba Va Gaa Va Qab (23b') 
as summation over discrete elements ready for use in a calculation then reads: 

N 

G{f,r') = g{r,r') + ^gif.rj) V{fj) g{fj,r') 

3=1 

N N 

+ T.T.ai^^''^^ ^(^"'^■) ^(^"':''^"'*) ^(^"'*) 9ir^,r') (23b") 
j=i 1=1 

Of course, summations run over all object cells here. No numerical advantage can be drawn 
out of r' = r in (23b"). With the same effort of making Gaa it can be used to evaluate 
maps of G{f, r) as well as plots of G(r, r') with r' fixed or even some function of r. 

II. 7. Conclusions and Outlook for section II 

We have presented a method to solve the problem of scattering of electromagnetic waves 
off an arbitrary distribution of dielectric objects, that is the exact evaluation of the field, 
especially in the near zone where higher-order multiple reflections can become important. 
Besides the field distribution we have obtained the Green's tensor characterizing the system 
independently from the form of the excitation. It represents the response function and also 
the density of states for supported electric fields. 

For a methodical introduction we have restricted our considerations to dielectric materials 
and the electric field. Without magnetic susceptibilities the magnetic field distribution can 
be calculated once having the electric field inside the objets by a formula like (9) with the 
magnetic background field and replacing g hy a. tensor including the conversion from the 
electric to the magnetic field by taking the rotation [13]. It is further possible to treat 
non-uniform magnetic permeabilities and even mixed systems with dielectric and magnetic 
objects [30]. The electric Green's tensor presented above is then paired by a magnetic 
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counterpart and genuine mixed response functions also exist. Whereas the calculation of 
the field distributions even for mixed systems is quite straight forward, the construction of 
the Green's tensor is more involved. It lives of the idea of handling one kind of objects first 
and then considering this setup as the background to include the other kind. There is no 
approximation or ranking in importance in this procedure. 

The discussion here has only considered finite objects in a homogeneous background 
as well as cartesian coordinates where vector and tensor components have been written 
out. Cylindrical and spherical coordinates are also commonly used [321 HI] and the Green's 
functins formalism has been developped for layered media f2E\ 115] . Besides wave- 
guide applications the use for modelling typical near-field optics experiments, where the 
microstructures to investigate are prepared on a substrate surface, lies in putting the influ- 
ence of this surface into a background Green's tensor [331 113] which is then implied the way 
we used g here. 

Details of applications of the Green's functions technique in electrodynamics to more 
complicated situations as well as beautiful results of corresponding experiments can be found 
in the given references. This text focussed on calculation techniques and further intended 
to give an overview of slightly different formal ways to calculate Green's tensors and fields 
of which either may be optimal for a specific problem. 

Appendix A: an unusual matrix inversion 

Suppose a complex quadratical matrix to invert is already given in the form M = 1 — M = 
1 — (m)jj or if it is not, we rewrite it like that. There is no restriction on the values of the 
numbers rriij. To get the inverted matrix proceed as follows: Of matrix M one by one take 
the diagonal elements ma and to all elements niab add mai{l — mi.;)~^mn,. After having 
worked through the matrix for one such ma, the changed matrix values have to be taken 
to do so for the next, also already changed, diagonal element. Obviously such steps are 
required for an x A^-matrix. This will yield (1 — M)^^ — 1, such that in the end 1 has to 
be added to all diagonal elements in order to obtain (1 — M)~^. For clearness we write out 
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the first two transformation steps of the matrix: 



mil "7,22 

m2i m22 



V 



mil + mii(l - mil) ^rnu mi2 + mii(l - mn) ^mi2 
m2i + m2i(l - mii)"^mii m22 + ^21(1 - mii)"^mi2 



^ mil + mii(l - mil) ^"^11+ mi2 + mii(l - mn) ^mi2+ ^ 



(mi2 + mii(l -mil) ^"^12)- 



(mi2 + mii(l - mil) ^"^12) 



(1 - m22 - '"T'2l(l - ?7T.ii) ^mi2) ^ (1 - m22 - ?7l21 (1 " ?7T,ii) ^mi2) ^ 

•(m2i + m2i(l - mii)"^mii) •(m22 + m2i(l - mii)"^mi2) 

m2i + m2i(l - mii)~-^mii+ m22 + m2i(l - mii)~-^mi2+ 

(m22 + m2i(l - mii)"^mi2)- (m22 + "^2i(l - muy^niu)- 

(1 - m22 - rn2i{l - mii)"^mi2)"^ (1 - m22 - "^2i(l - rnii)~^mi2y^ 
\ -(^21 + m2i(l - mii)"^mii) -(^22 + ^21(1 - mii)"^mi2) 



(24) 



The inverted matrix (1 — M) ^ can be represented as a geometric series: 

(1 - M)-^ = 1 + M + + . . . 



(25) 



Truncating and using the sum from the right side is only possible if the series converges 
whereas the closed form on the left is valid in any case. In contrast to the infinite sum on 
the right side of (25), our inversion procedure consists in a finite number of steps of adding 
contributions to the matrix elements. Nevertheless, (25) tells us that the invers (1 — M)~^ 
is the sum of all powers of M and thus each element in row a and column b must be the sum 
of all possible products ma„jm„^„2"^n2n3 ■ ■ ■mn^ mjnT.njb with any number j of inner indices 
including none. (Diagonal elements get an extra +1.) There are N different indices rii and 
they may repeat, of course. Considering that (1 — mjj)"^ can also be written as 



(1 - ma) ^ = 1 + mjj + m-j + 



l + mii + mumu + 



we see that the first step in (24) adds to each matrix element ruab the sum of all products 
maimiimii . . . mumu,. In these at least one pair of indices 1 is squeezed between a and 
6 as in maimn,, the contribution niab was already there. In the second step all products 
mav-^m^^t^^ . ..m„^2m2vf3 ■ ■ ■ my^_^y^mi,^b with V1V2 . . . I'j-ii'j every possible sequence of Is and 
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2s will be added. The products with only indices 1 between a and b were there before. In 
the third step every sequence of indices 1, 2 and 3 with at least one 3-hnk is added. And so 
on until in the end at each matrix position between outer indices a and b we have created 
all possible sequences of an endless game of dominos rriij with numbers i and j from 1 to 
A'". This argument was to proove that the result of (24) indeed gives (1 — M)~^ — 1. We 
calculate a finite number of (1 — mu)'^ or products mai{l — mii)~^mn). The sequence on 
the right side of (25) need not converge and the original entries in M need not at all be 
small compared to 1 in their absolute values. There is no approximation in the sense of a 
perturbation theory. If the matrix 1 — M is degenerate, the failure of the inversion will be 
noticed when a value 1 — ma becomes zero at some step. Not to confuse notation, remark 
that in (24) and in products in the text like maimn . . . rriib letters m meant the original 
matrix entries whereas in expressions 1 — ma, (1 — ma)'^ and mai{l — mii)~^mn) we referred 
to the entries at the respective step of the matrix transformation. 

In the application from the main text i — 1,...,N enumerates the object cells. At 
position in (1 — gV)~^ when expanded into a series having every possible sequence 
5'(^ni, ^«2)^(^n2)5'(^«2! ^n3)^(ni3)fi' • • • 9(^,^-1, rn,)V {frij) shows that the resulting field at any 
place (inside or outside the objects) is the interference of the background field and the fields 
reradiated by all the object dipoles having undergone every possible scattering path between 
the objects (Fig.l). A complication in the electrodynamics application at this stage is the 
fact that each matrix element m actually in itself is a 3x3 matrix indicating the effect of 
three field components at one place onto three field components at a another place, g on the 
discrete space of the objects can be written blockwise with a scheme of quadruple indices 
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{fi,x;f2,x} {fi,x;f2,y} {fi,x;f2,z} 



{fi,y;fi,x}({ 



rL,y;ri,y}) {fi, y; n, z} 



{fi,y;f2,x} {fi,y;f2,y} {fi,y;f2, z} 



{fi.z-.fi. x} {fi ,z\fi, y} (^fi.z:fi,z}~^ 



{fi,z;f2,x} {fi.z;f2,y} {fi, z;f2, z} 



i_ 



J 



{f2,x;fi,x} {f2,x;fi,y} {f2,x;fi,z} 



{f2,x;f2,x 



T){f2,x;f2,y} {r2,x;r2,^} 



{f2,y;fi,x} {f2,y;fi,y} {f2,y;fi,z} 




{r2,z;fi,x} {r2,z;fi,y} {r2,z;ri,z} 




V{fk) in the product 5'(fj, ffc)y(ffe) just multiplies the respective column. 

One could use the inversion procedure working off the diagonal elements marked by ovals, 
requiring 3N steps then. However, the process is equally applicable to 3x3 blocks as marked 
by the dashed rectangles, since by its dimension the whole matrix can be divided up into 
3x3-blocks. Then (1 — mu)'^ means the inversion of a 3x3 matrix and mai{l — 'mii)^^mih 
to update the blocks means the product of three 3x3 matrices. These operations should be 
programmed as elementary procedures. 

The given procedure to invert a matrix can become of advantage if for sparse matrices 
conventional routines run into numerical difficulties because of many zero values. Besides 
that, it can be adapted to become quite efficient if only parts of the inverted matrix are 
needed or for symmetry reasons it is known that blocks or patterns of matrix elements 
vanish and will stay zero throughout the inversion. Although for the calculation of the 
Green's tensor a slightly modified procedure is applied that directly optimizes the numerical 
solution of the Dyson equation (see appendix B), the matrix inversion was discussed here, 
because it may be used for more general purposes and in other contexts as well. 
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Appendix B: Calculating the Green's tensor 

In the following instructions are given how to calculate the Green's tensor 

G = (1 - gVy^g = g + gVg + gVgVg + ... (21) 

being efficient in the way that finally only G(r, f) of equal arguments on the mesh points 
of map B need to have the correct values [501 1^ . For a start nevertheless consider that 
G(rj, fk) with arguments fi and from the big set of all object and map cell midpoints will 
have to be the sum of all products 

)V{rnJg{rm,rn^)V{fn^)g(f ) • . .V{rn,.r,rn^)girnj,rk). 



Tux-i 'f'n2i ■ ■ - I "f^rij cau ouly bc object cells and any sequence of them has to be created, including 







AB 


the empty one with no V giving the term g{fi,rk)- Initiate a matrix 




BB 




^BA 


BB 



- for simplicity call it G from the beginning - with line and column arguments running 
over all object cells in A and all map cells in B with g{ri,r]t) in each 3x3-subblock. In the 
5-B-quadrant only diagonal blocks g{ri,fi) will be needed, however. Work off the diagonal 
subblocks through the AA-quadrant. For the nth one prepare the inverted 3x3-matrix 
(1 — G{fn,fn)V{fn))~^ and to all entries from line n + 1 and column n + 1 on add the 
product given below. 

G{n,fk) ^ G(f„ffe) + G{n,fn) (1 - G{n^,f^)V{f^))-' G(f„,ffc) (26) 

In the following step use the updated entries in the above recipe. In the nth step you need 
not update entries in line n or above or in column n or to the left of it, because these will not 
be needed as multiplication factors (^(rj, rm) and G{fm, r^) for m > n any more. The parts 
of quadrants AB and BA remaining after these restrictions have to be changed by (26) for 
every n. Merely diagonal 3x3 blocks have to be done in the SS- quadrant. The first step 
adds all products consisting of any non-zero number of factors V{ri) between gs, the second 
step all products eventually containing V{fi) and one up to any number of V^(r2), and so 
on. The procedure is finished after fvi has run down the diagonal of the AA-quadrant. All 
sequences of multiple scattering from the objects are then included. 
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We could have filled the whole SS-quadrant with initial qbb 3x3-blocks and updated 
the entire SS-quadrant in each step. Like in the matrix inversion procedure from appendix 
A we should further have worked down the complete diagonal and done (26) for every fVi 
from region B as well. There will, however, be no additions as V{fn) = for fn in B 
(even if a map cell accidently coincides with an object cell). This argument also reveals 
why non-diagonal elements in BB do not have to be evaluated. They can never appear as 
multiplication factors G{fi,fn) or G{fn,fk)- The BB diagonal is the resulting (^(f, r)-map 
we wanted. 

Initiating the matrix by g and introducing the V on treating the respective diagonal 

element makes the outcome of the procedure directly (1 — gV)~^g compared to initiating 
the matrix with gV and as an intermediate step obtaining (1 — gV)~^. 
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III. ALL-ORDER QUANTUM TRANSPORT 



It is demonstrated how the transport problem for two open free-electron gas reservoirs 
with arbitrary coupling can be solved by finding the system's Green's function. In this sense 
the article is an introduction on Green's functions for treating interaction. A very detailed 
discussion of the current formula is given on an elementary basis. Despite formal resem- 
blances the stationary transport situation, however, differs in its nature from introducing 
coupling between energy levels in a closed system where then the interest lies in modified 
eigenvalues and eigenstates. 

III.l. Introduction 

By modern lithography techniques so-called point contacts [17] can be arranged between 
conductors. These have sufficiently small dimensions such that electronic modes get quan- 
tized. However, coupling across such constrictions need not be so weak as to be described by 
a small tunnel probability, but can be influenced by coherent interference of multiple reflec- 
tions. Point contacts can be obtained by indenting STM-tips into some material |18], elec- 
tromigration |49j or the break-junction technique [50]. Whereas for constrictions imposed by 
gate electrodes to a two-dimensional electron gas in semiconductors [51J one observes quan- 
tized conductance values in the sense of fully transmitting or totally switched-off modes, the 
application in mind behind this work is the type of connection like the single-atom contact, 
characterized by an ensemble of channels [52] , which can also have intermediate transmission 
amplitudes between zero and one [53]. Even with some fully transmitting modes the contact 
bears a resistance in the order of the quantum resistances Rk = h/e^ = 26kQ [31], such 
that viewing the system as a left and a right side with some interaction is an appropriate 
picture. Furthermore contributions from several channels just add in the current. Although 
a transmitting channel in a point contact is the application in mind, this article shows how 
to set up a general procedure to solve the problem of transport for two open reservoirs with 
more or less strong coupling between them, and thus how the Green's functions formalism 
from section I is implied in an area of current research [HSl ES] • 
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III.2. Green's functions formalism 



(1) 



As a preparation, consider systems like, for example, the bulk material on the left or the 
right side (Fig. la) without coupling. For these we suppose that we know the Hamiltonian 
and the wave function at any given energy hu satisfying the Schrodinger equation. 
Denoting the Hamiltonians Hll = Hrr = (doubling the index makes sense later), the 
wave functions — — ^"^^ the time derivative as dr, the Schrodinger equations for 
both sides - for the moment just formally put into matrix form - are 

\hdr -Hll \l ^0 (r)\ _ 1 

Corresponding to this differential equation we have the Green's functions equation 

ihdr -Hll \ IgLLir, r') ^ \ ( '^(^ " 

ihdr-HRR) \ 9rr{t,t')) \ 5{t-t')^ 

(2) 

If we knew g — Qll — Qkr, we could immediately also give a solution to (1) with a source 
term added 

Uhdr-HLL \ f^lir) 







ihd-r — H 



RR, 




(3) 



namely 



dr' 









(QLir'n 




[QRir') 



(4) 



9LL[r,T ) 

gRR{r,T' 

However, we are more interested in the solution when an interaction between left and right 
is present, expressed through coupling Hamiltonians Hlr and Hrl, 




(5) 



ihdr - Hll -Hlr 
—Hrl ihdr — Hrr 

and we shall refer to this case by ^ without upper index. Introducing the coupling, of 
course, was the motivation for writing the Schrodinger equation in matrix form over the site 
space consisting of L(left) and R(right). Putting the coupling terms on the right side of the 
equation, they mimic a source 



ihd-r — H, 



LL 





ihdr - Hrr 





(6) 



31 



and following (4) the solution can formally be written as 

.9LL{r,T') \ I Hlr 
gRR{r,T') [hrl 





(7) 



obtaining the implicit Lippmann-Schwinger equation for the wave functions. Although for 

the coupled system we do not expect the eigenvectors of 1 to be one with only 

an upper component localized on the left and one with only a lower component localized 

on the right, it is convenient to denote these solutions as vectors , keeping indices 

L and R. In contrast to section II, where the Lippmann-Schwinger equation was solved in 
a discretized form to obtain field distributions, for our purposes here (7) will merely serve 
as a formal step in the derivation of the Green's function. What we are precisely interested 
in from a physical point of view is the current that will flow between the left and the right 
side, and not necessarily explicitly evaluating the wave functions. 

The Green's function G of the coupled system is inferred from an equation analogous to 
(5), namely 

Uhdr - Hll -Hlr \ I Gll{t, t') Glr{t, t')\ ^ I 5{t - r') 
-Hrl ihdr-HRRj \Grl{t,t') Grr{t,t')) \ 5{T-r')^ 

(8) 

G like the Hamiltonian is a full matrix in site space. With G, a solution of 



-Hrl ihdr-HRR) U^^(r) 




(9) 



would read 



(10) 



^ 1^^, (GllM GLR{ry)\ /g^(r') 
,^^VV ' ^ \Grl{t,t') GRRiry)) \Qr{t') 
where the index cq stands for coupling and source. As explained in section I, we could also 
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have set up a Lippmann-Schwinger equation for as 

V/(r)^ Z"^^, fgLL{r,T') \ (Q,{t') ^ 



dr' 



/ 



(11) 



Now from inserting (10) into (11) and for Q choosing some 5(t' — tq) either in the L- or in 
the R-component the Dyson equation for G is obtained: 



+ 



Gll(t, t') GUr. ] ^ / ^ll(t, t') 

G«l(t,t') Grk(t,t')/ \ ^rr(t,t')^ 

'^?LL(r,r") Wo Ig^t^t') G^^" ,t') 

5ilR(T,T")j \iinL I \Grl{t",t') Grr{t",t') 



dr" 



(12) 



III.3. Explicit Green's functions in time and frequency domain 

The reservoirs on the left and right being bulk metal, the differential operator of the 
homogeneous differential equation (1) is given by the free-partical Hamiltonian 

ihdr - //° = ihdr + — A (13) 
2m 

and the corresponding wave functions are — Q±ikr-uvT according to the dispersion relation 
E — hjj — ^k'^- Looking for a Green's function at given frequency u, in (2) we can replace 
by fko and thus have to solve 

{ihdr - hj) 9{r, r', u) = S{t - t') (14) 

We need an expression, the derivative of which produces 5(r — r'). One easily verifies that 
(14) is satisfied by either 

/(r, r', a;) = ^(r - r') e'^'^^^-^') (15a) 

g^ir, r', c^) = ^ 0{t' - r) e'"^^^-^'^ (15b) 
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FIG. 4: A left and a right bulk reservoir, uncoupled (al) and coupled (a2). (b) Contact as embedded 
in a circuit, (c) Fermi levels left and right with applied voltage. The dotted line marks a fixed 
energy over both sides. 

The retarded function only exists for t > t' and the advanced function g"- for t < t'. In 
(14) for simphcity we restricted ourselves to a single energy and therefore (15) still contain 
a; as a parameter. For complete Green's functions in the time domain these terms have to 
be multiplied by the density of states ^{00) and integrated over energy. Later in transport 
we shall be interested in the amount of charge transferred, not resolving any more, which 
energy levels contributions came from. Thinking physically of a small contact between two 
metallic leads one might argue that in the constriction transverse A;-vectors are quantized 
and a one-dimensional continuous density of states remains. Nevertheless, for not too high 
voltages only a certain energy range around the Fermi energy will play a role in transport 
and thus setting the density of states equal to its value at the Fermi energy is a good 
approximation. Anyway, a constant density of states V^u) = V with occupied states below 
and empty ones above the Fermi level for each side, left and right, shall just enter our model 
here as an assumption. We are led to the following representations of the retarded and 
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advanced Green's functions: 

= ^ j duj g^to) e-'"("-"') with g'-iu) = -i (16a) 

g''iT,T') = j duo g\uo) e-'"^^^-^'^ with g\uj) = i (16b) 

Taking out the factor "D/Zi, we get the dimensionless functions g'^{oj) = —i and g°'{uj) = i in 
frequency space. Remark that their deduction here did not consist in calculating a Fourier 
transformation. The phase factor e~'^^^'^~'^'^ was already there in (15). The tu-integrals from 
(16) will be implemented as a useful representation of g"^^"" still in the time domain. The 
6'-functions have been deliberately skipped after the =-signs. We shall however see that 
g^ and g°' finally only appear with the correct relation between their first and second time 
argument. The a;-integral should indeed be understood as summing over all energies and, 
even if g^^"'{uj) is a constant, on no account be interpreted as this constant times 2ii5{t — r'). 
g-ia;(r-r') ^ ^-iujT _ (^^-iujr' y ^ |vl/0(r) >< ^°(r')| Is the projectiou-like couvcrsiou of the phase 
unit vector of an oscillation from time r' to time r. Left(L) and right(R) are distinguished as 
the origins of wave functions constituting our basis of states, but the junction is considered 
point-like, that is r* = for both left and right and thus e^^^^ = 1 drops out in an alike 
combination of \E''^s even from L and R. The two-fold time argument suggests to take (16) 
formally even as double Fourier transform (however with different signs in the exponentials 
with r and r') 

^ //'^(r,r') = jdu^ j duj2 g'/%oo,,U2) e''^^^ e'^^^' (17) 

with g'^^°'{uji,uj2) = =F^5(co'i — 102), however, because different energies stay independent. 
(More or less guessing the Green's functions g was easy in our normal-conducting simplest 
model here. Generally, it has to be found as the solution to an equation like (14) with the 
differential operator from the uncoupled system's Schrodinger equation and an elementary 
perturbation 6. In the superconducting state, for example, the two-particle interaction in 
the Hamiltonian forces the Green's function to include Andreev reflection |5T], and the result 
is such that g{uj) does not just vanish in the gap of the quasi-particle density of states.) 

Before we can solve (12) for the coupled system's Green's function G, we have to specify 
the coupling parts of the Hamiltonian Hlr and Hm. Like with the density of states the 
simplest model will assume that the coupling is energy-independent and described by a 
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constant (real) interaction energy W. The reservoirs are unaltered by transport between 
them. Due to the apphed voltage (Fig. lb) incoming charge carriers (from the left on the 
right) are led away and outgoing ones (on the left to the right) get replaced. In the contact 
wc do not allow relaxation or other energy-changing processes. Wc choose the repective 
Fermi levels as zeros of energy on either side (Fig.lc). An electron going from right to left 
has to strip off its phase e"*^^"^ and aquire e~^'^^'^ to fit in on the left. An analogous argument 
holds for transitions from left to right, and therefore the coupling terms are 

Hlr = (JLRir) = W e-^('^^-'^«)^ = W e'"^^'^ (18a) 
Hrl = (Trl{t) = W e-^(^«-'^^)^ = W e-'^^^'^ (18b) 

The phase factor gives a time-dependence to Hlr and Hrl, but ujr — ujl — eV/h is 
indeed independent of energy. For G we make an ansatz hke (17) as a two-fold Fourier 
representation: 

^ G{t,t') ^ j dwi j du2 G{ui,uj2) e-^'^i^ e'""'^' (19) 

From (12), which is valid for either advanced or retarded functions, as an example, we pick 
the upper left component of the 2x2 matrix in LR-space and insert (17), (18) and (19): 
(Multiple integral signs are skipped from now on.) 

Gll{t, t') ^ gLL{T, t') + j dr" gLL{T, r") Hlr{t") Grl{t", t') ^ 

j dooLi dujL2 Gll{ojli.ojl2) e"^'^"^ e^'^""' = 

j diVLl dUL2 gLL{i^LUi^L2) e'"'^'^ c"^^^^' + 

dr" djjjLi douLs douR du}L2 5'll(<^li, c^ls) e"''^-^!'" e^'^^^-r" . 

• W e^^^^"/'^ ^ Grl{u;r, ujl2) e-'"«"" e^'^-^' (20) 

The integral over r" produces 2Tr5{u)L3 + ^ ~ '^r) in giutLiji-^Ls) there is 5{u!li — ^^ls) 
anyway, such that the last term of (20) becomes 



V 

~h 
V 



— / dULl duJL2 g{ojLi) t Grl{ojli + -^,0Jl2) e 



with t — WV/h. Strictly speaking, if (20) were for the retarded function, in the last term 
there would be 9{t — t") from g'£^ and 9{t" — t') from G^^, and if it were for the advanced 
function, 9{t" — t) and 9{t' — t"), such that the integral over t" only exists between r and r' 
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instead of having minus and plus infinity as fimits. However, with culs and cur running over 
any value, one can argue that even with a finite T"-integral the only remaining contribution 
stems from cui^3 + ^ — cuji = 0. A discussion of time ordering will again appear in section 5. 
V is the density of states per frequency interval, dividing by h makes it the number of states 
per energy interval. W is an energy, t can be understood as a dimensionless transmission 
amplitude. Now we set up the convention that all frequency arguments of g and G are 
written with respect to the left zero level and for the case that they correspond to the right, 
that is an R-index, it is understood that eV/h is added. (20) has to hold for any r and r' 
and from comparing Fourier coefficients we get 

Gll(<^li, ujl2) = gLL{uJLi) S{ujli " <^L2) + 5'll(<^li) t G rl{uj LI , L2) (20a) 

We could have inserted Grl = QrrHrlGll and again replaced Gll = Qll + Qll^lrGrl 
and so on. Instead of an implicit equation for G this would have led to an infinite series (see 
section I): 

00 / CO \ 

G = gJ^^'^aT = ] 9-9 + 9(^9 + 9(^9(^9 + ■■■ (21) 

n=0 \n=0 / 

(21) is written for whole matrices in LR-space, calculations like (20) can be done analogously 
for Glr, Grl and Grr. (21) can be read as equation in the time domain. Then a is 
the matrix consisting of Hlr and Hrl and each multiplication of two following gs with a 
in between means an integration over time. However, (21) is as well valid as relation in 
frequency space. In this case a is the 2x2 matrix with just t as off-diagonal elements. Like 
we have seen through evaluating the r"-integral in (20), in (21) each connection Hlr/rl 
from a passes the frequency argument from the g in front to the g behind. And as g can 
only have two identical frequency arguments, no uj different from the first can ever appear, 
such that ujli = (^l2 — ^ and G also effectively is a function of only one frequency argument: 
G{uji,uj2) — 5{ijJi — uj2)G{uJi). With Green's functions of a single frequency argument (20a) 
and its analogues for the other three components in LR-space become a simple algebraic 
equation: 

'Gll{uj) Glr{uj)\ ^ I gLiH \ / gLLH A ^llH Glr{u;) 

Grl{u;) Grr{uj)) \ gRR{uj)) \ gRR{uj)) \t \Grl{uj) Grr{uj) 

(22) 

(The fact that even the Green's function of the coupled system turns out to be a function of a 
single frequency argument is a special feature of our simple model for the normal conducting 
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case. In the extension of this model to superconducting reservoirs and transmission processes 
including Andreev reflection, G becomes a function of two frequency arguments, the second, 
however, restricted to values differing from the first by an integer multiple of eV/h, such 
that effectively there is one continuous and one discrete frequency parameter [58j.) Here, 
with g'^'/^-^uj) = =Fz, (22) is easily solved and G{uj) comes out independent of frequency, too: 

r/a 

Gll Glb. \ , . 

Grl Grr I 

1 / -To _+ \ 

(23) 

The same result could have been obtained from (21) by writing out a few more of the matrix 
multiplications and using the formula for the geometric series in each element. We shall need 
two further types of Green's functions. T will be introduced in the next section and 
and G'^" when calculating the current. 

III. 4. Transfer Green's functions 

In (21) there was a sum of products of arbitrary many factors g and a with outer factors g. 
"Product", of course, except with the Green's functions taken of a single frequency parame- 
ter, in the time domain or with two-fold frequency dependence still meant a convolution-type 
integration over inner arguments. In analogy we define the sum of products with outer fac- 

10 t 

tors a (as integrals in the time domain or just algebraically with g{uj) and a 

\t 

oo / oo \ 

T = a + aga + agaga + ... = a Y,i9< = J2^cTgr a (24) 

n=0 \n=0 / 

Whereas all contributions to G in (21) began and ended with staying some time in a reservoir, 
described by (yf - at L from r to r" for a start in the last term of (20), for example - each 
term of T in (24) begins and ends with a transition a and further contains at least one such 
hopping across the junction (which g does not). Therefore we call T the transfer Green's 
function. The same as for G, the relation between T as a function of times and as a function 
of frequency is given by 

I Tjk{t, ^') = ^J duJ TjKiu) e-'^^^ e'^^^ (25) 
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where cuj — cu ii J—L and cuj — cu + eV/h if J=R and the same for ujk- (Taking out 27r of 
V in (16) was a convention. The prefactor of the a;-integral for T follows from consistency 
requirements. (19) for G was in complete analogy to (16) for g. T with two time arguments, 
however, has a little different character from a with just one time parameter.) Alternatively 
to (24) T could be defined through its link to G 



aG = Tg or Ga = gT 



(26) 



Be careful that replacing one by the other can introduce another internal time integration 
as, for example, a{r) G{t,t') = J dr" T{t,t") gir'^r'). (24) and (26) hold for retarded 
and advanced functions. From (24) it is immediately seen that T like G satisfies a Dyson 
equation 

G^g + gaG and T ^ a + agT (27) 
and even the complementary forms 



G^g + Gag and T = cr + Tga 



(28) 



are analogues. In Fourier space, hke (22) the T-equation (27) is an algebraic equation and 
the solution like 



G^il-ga)-'g is T ^ {1 - ag)-^a. 



(29) 



Inserting g and a explicitly for our model we get 

r/a 





especially 



r/a 



-irja 



t-t^ + t^ - ... 



(30) 



(31) 



LR — -^RL I J^f 'i 

Whereas t is the single hopping amplitude, Tmij^ is a renormalized transfer amplitude. 
One may wonder why a model for transport could not have been set up adding amplitudes 
for transfer processes of all orders, as the interaction (18) seems to be introduced the way it 
is just in order to result in powers of t. However, g'^/"'{uj) — deduced from the Schrodinger 
equation is decisive for the signs in (30) and (31). One may wonder that multiple reflections 



39 




FIG. 5: Transfer processes from left to right of different order. 

are not added as t + + + . . . = P^^ocesses of different order (Fig. 2) are not 

independent, but interfere. Trl is the transfer amphtude per single electron supplied on the 
left by the voltage source. But for every electron that goes over to the right in an nth order 
process {n odd) with weight there is one that has hopped once more to the right and back 
(n + 2 order process) and thus is not to be newly supplied, but to be again sent through 
the junction. The amplitude t is renormalized by 1 + as denominator. However, such 
interpretations of quantum mechanical amplitudes are precarious, and the full conversion of 
t to a transmission probability will be established later. 

It is quite instructive to solve (27) in a slightly different way than done in (30). Firstly, 
for the four components in LR-space we have 

Tll = (^lrQrrTrl Tlr = ctlr + c^lrOrrTrr 
Trl = CFRL + (TrlQllTll Trr = (TrlOllTlr (27a) 
Inserting these into each other, for example, an equation for Tlr alone is obtained: 

Tlr = <JlR + (^LRgRR<^RLgLLTLR (32) 

This implicit equation is the basis for calulating the transfer Green's function in more com- 
plicated cases than discussed here [5HI EH], like for example the superconducting junction. 
In our model, inserting aLR = ctrl = t and g^/l{uj) = g^^^{ijj) = into (32) immediately 
also leads to Tl^^{uj) = 
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FIG. 6: Illustration of the Dyson equation for the transfer Green's function. 

More easily than the Dyson equation for the ordinary Green's function G, the one for 
the transfer Green's function T is illustrated as is done for the LR-component in Fig. 3. 
(Normally indices are read from right to left such that T^r is considered a transition from 
right to left, but it does not really matter whether they are interpreted the other way round 
as in Fig. 3. The actual sequence of what is earlier or later in time will be discussed when 
calculating the current in the next section.) Fig. 3 demonstrates the implicitness of the Dyson 
equation: Any transition from left to right is either a single transfer or an electron hopping 
to the right and back followed by any process beginning on the left and ending on the right, 
no matter what happens in between. This last part by definition is the sane as the other 
side of the equation, namely Tlr. 

III. 5. Calculating the current 

From the Heisenberg picture of quantum mechanics we know that the time derivative of a 
not explicitly time-dependent operator A is given by the commutator with the Hamiltonian 

iA=UH,A] (33) 
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The operator of interest here is the projector on either side of the junction 



PL 













2 




or PR 



in 


><n 








2 




(34) 



As explained earher, with the junction couphng left and right together, the solution 
is not limited to one side, however, the projectors take out the respective part: 

and 








(34a) 



Pl and Pr are proportional to the amount of charge on the left and on the right side. Their 
time derivatives have equal absolute values, but opposite sign and represent the current. 



^ dpL ^ ^ dpR 



(35) 



dr dr 

(< I and I > are used for bra- and ket-states. Here < > means the expectation value, 
of course.) e is the charge of an electron and the sign of I can be defined arbitrarily. We 
choose to do the calculation with p^ and evaluate the commutator with the Hamiltonian: 



[H, Pl] 



Hll Hlr 

Hrl Hrr 





(36) 



Putting together (33) and (35) we obviously need the expectation value of the operator 
[H,pl\. [H,pl\ shall be called (Xc- Even if (33) stems from the Heisenberg picture, it is 
written in such a way, that the right hand side is to be evaluated in the Schrodinger system 
with time dependent states, and we shall here change to the interaction picture [Mj for the 
calculation. With standing for complex conjugation as well as transposition from column 



to line vector, the value of 



dpL 
dt 



in state 





1 M 








at time r is given by 



< ^°(ro)|T exp(- J i a^(r') dr') i[H,pL]n T exp(- 



-i)an{r') rfr')|^°(ro) > 



(37) 
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Hlr{j) and Hrl{t) denote the time occurence of r from (18) which has to be considered as 
still belonging to the Schrodinger picture in our case here. T means time ordering [HI] and 
T anti-time ordering. in the uncoupled system, of course, also stands for a two-vector 

with left and right component. Replacing by '^^{tq) in (37) we took out both the 

coupling as well as the time dependence from the states. Accordingly the a-^ in the integrals, 
in contrast to the Schrodinger-picture a used in the preceding sections, in a Heisenberg way 
have to include the time dependence of the uncoupled states. o"(r) could already be written 
as 

^ i^L/?(r)\ fe-'"^^^ \ M/\ /e*^^^ 
Hrl{t) J y e-'^R^j \ w Q j y Q e*^« 
or cr(r) = |\1''^(t) >W< ^°(t)| for short, where bra, ket and W still mean the respective 
matrices. However, this decomposition might rather be confusing and will not be used, 
anyway. The translation to the interaction picture is the following: 




\ /e"^^^ 








(38) 



or (T7^(r) = |\E'°(ro) >< ^"(t)! cr(r) |^°(t) >< ^(to)]. The Heisenberg picture always refers 
the operator back to the undeveloped state at tq. 

Analogously will have to be extended by the time-dependence of the uncoupled states. 
The meaning of the time-integral over a-^ as an exponential is best explained by writing 
explicitly: 

^ \ 1 /"^ 1 1 f^n-l 

T exp(- 1^ (->) drO = E fi X ■ ■ R ■ 

{-i)(rn{Ti) . . . {-i)an{Tj) . . . (-Of^w(^n) (39) 

where all arguments ti, . . . ,Tj, . . . ,Tn have to lie between tq and r and all products from 
none to arbitrarily many factors a-^ have to be added. With time ordering, furthermore 
Ti > . . . > Tj > . . . > r„ is imposed. Without T the scheme for the exponential series 
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FIG. 7: The Keldysh contour illustrating the development out of the uncoupled system's states of 
the bra < *(r)| on the minus and the ket |*(t) > on the plus branch. This view picks out the 
single transition at time r and visualizes the calculation of < *(r)|(7(r)|*(r) >. Each point on 
the contour represents an LR- or RL-transition described by a, each line segment the development 
of the wave-function phase which is given g. 

e"" = Er=o ^ in (39) would produce J2n^A Iro dn ■ ■ ■ \ f^^ drj . . .\ f^Jr^ . . . . But then 
for every fixed set {ri, . . . , Tj, . . . , r„} there are n\ permutations for the time values to occur 
such that time ordering and restricting the arguments to mutually exclusive intervals as 
done in (39) cancels the factorial denominators. Writing out the anti-time ordered part in 
the same way as shown for the time-ordered part, the expression from (37) becomes 

< *«(to)| *°(to) >'-< ^\Tm)\a{Tm)\^\Tm) >< *°(to)| ... 

. . . |*°(to) >'-< ^°(r,)|a(r,)|^°(r,) >< ^°(ro)| ^^r,) > ^ < 

a{n-l)\ . . . W{na)mna) >< *°(to)| *°(to) > ^ < *"(T)|<7e(T)| 

*°(r) X *°(to)| *°(ro) > (^) < ^'{nrMn,)\ . . . 

...|<7(r,_,)|*^«(r,_,) X ^"{to) |*Vo) > (y) < *V.)k(r,)|*°(r,) > 
<*°(To)|...|a(r„)|vl/0(r„)><*°(ro) |*°(ro) > (37a) 

where Tq < < . . . < < . . . ria < r, Tq < r„ < . . . < r^- < . . . < Ti^ < r and integration 
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over all time arguments except r and tq is understood. Like in (39) we mean the sum of all 
products with arbitrary many inner development factors like of and of Tj. Spaces 
and different font types are just used in (37a) to recognize sequential cr-^-parts like (38). 
Inner development factors like \['° of and Tj as well as '^^{t) and the inner ^°(ro) have 
to sum over space or the basis of states and therefore should be thought of as matrices as 
in (38). We have not yet decided which components outer bra and ket with *°(to) will 
project out. All < ^°(to)|^'°(to) > are the unit matrix and drop out (^° being given by a 
non-normahzable e-function does not pose a problem here). Then |^°(Tjk) > | < ^°(Tfe_i)| 
with Tfc < Tfc_i is recognized as g°'{rk,Tk-i) and |\E'°(rj_i) > < ^'°(t,)| with r,_i > tj as 
g^(Tj^i,Tj). These are 2x2 diagonal matrices. With (21) in (37a) the whole sequence 
from |\l'°(rm) > to < ^°(t)| can be replaced by G"(rm,r) and the long part from |^°(t) > 
up to < ^°(r„)| is just C(r, r„). The complete operator between the outermost < *°(to)| 
and |*°(to) > in (37a) therefore becomes 

|*0(to) >'-< ^°(r^)| [1 S{Tm - r) + a{Trn)G''{Tr^,T)] a,(r) 

[lS{rn-r) + G^{r,rn)a{rn)]\^\Tn) >< *°(to)| (37b) 

The 1-contributions stem from the cases where there are no factors with ^^{rk) or ^^(tj). 
We regard the coupled system as having developed out of the uncoupled system, but we 
are looking for a stationary state. To achieve this, the coupling has to have been turned 
on infinitely long ago, thus we let Tq — > — oo. The sequence of time arguments is usually 
represented on the so-called Keldysh contor (Fig.4). We still have to take the expectation 
value of (37) or (37a). This means summing over the basis of states for the outer ^"(tq). 






right for each frequency lo. We shall discuss the occupation or emptiness of states shortly 

after having worked off some further more formal points. Summing over 

the outer ^^(tq) will return the trace of the matrix M given in (37b). 

1 \ (mil mi2\ (l\ (o\ ( mil mi2 . , . ^ ,^ 

+ I I I = IV M = mil + m22 

/ \ m2i m22 1 \0} \1/ \ m2i m22 




— V 

M M 




Only being interested in the trace as a result, in a matrix multiplication the order of factor 
matrices can be changed cyclically. As in (37b) every bra and every ket as well as each 
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operator part between the two | is a matrix we can rotate factors to obtain 

a,(r) [1 5(r„ - r) + G"^(r, r„)a(r„)] |^°(r„) > • 
< *°(ro)|^°(ro) >'-< ^\t^)\ [1 (5(r,„ - r) + a(r„)G"(r„, r)] (37c) 

< \l>°(ro)|\I/°(ro) >= 1 again drops out. |\I/°(r„) > | < ^I/°(rm)| is of the structure g, 
however, with no restriction as to which argument r„ or is earher or later in time. We 
define this new type of Green's function as 

9^-ir,r',co) = le~^-^^-^'^ = { (40) 



(An eventually ill-defined single point r = r' is irrelevant for later integrations.) As the cases 
T > t' and T < r' are mutually exclusive the function can also be given by g^^{T,T') = 
g°'{T, r') — g'^^T, r') (conclusions on g^^{uj) from that are risky, to my opinion, though). We 
have thus deduced the current formula 

dTn / dTra [(l(5(r„ - t) + (t, r„)(T (?:„)] • 
■ ^+"(t-„, r„) [16{t^ -t) + a(r„)G"'(r^, r)] } (41) 

(like discussed in (20) it is rather irrelevant whether the upper integration limits are set as 
r or oo) or 

/(r) = -e Tr {a,{r) G'+-(r,r)} (41a) 

if is defined as 



G+-(r,r'l 



j dn j dT2 [l(5(ri - r) + G'\t, T^)a{n)] g^- [r^, r^) ■ 

[15(r2-r) + a(r2)G'"(r2,r')] (42) 



or 

G+- = (1 + GV) g+- (1 + aG") (42a) 

in short notation. The current only needs G^^ of two identical time arguments. We 
shall need the Fourier representation of g^~{ijj), but none such for G"*"". Although in our 
simple model for the normal-conducting contact the current will come out the same for any 
r, /(r) generally does depend on time (our expectation value is an ensemble, not a time 
average). For the superconducting case the ac parts [58] like the Josephson current are 
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included in the expression (41), although one is mostly interested in the contribution in I 
that is independent of r and gives the dc part. 

Using g~^~{Tn,Tm) without a third argument cu in (41) we understood integration over 
frequency like in (16). But this point requires more care as it has not been taken into 
account so far whether states are occupied or empty at tq — > — oo. Writing out the trace 
from (41a) in LR-components reveals four contributions 



It is clear that in contrast to the other appearing a, in ac from (36) there is an additional 
relativ minus sign between the LR- and the RL-component. Obviously the net current is 
the difference between the current from left to right and the current from right to left. In 
the first - as well as the second - term in (43) a transition from R to L at r, the time 
argument of a^, is picked out to be counted for the current (sec also Fig.4). With g^^ at R, 
the charge carrier is supplied from a state originally located on the right. It is an electron if 
the energy lies below the right-side Fermi level. Until r the state has evolved to again be on 
the right. The originally left behind empty state at the same energy below the Fermi level 
must have evolved to be on the left at r such that by the R^L transfer the electron can 
go into it (Fig. 5a). Or you might say that the plus and the minus branch of the Keldysh 
contor represent two possible parts in the evolution of an original wave function \E'5j, between 
which there is a non- vanishing matrix element of the operator ac — [H,pl]. However, still 
regarding the first line of (43) , there is the further possibihty that an empty state from above 
the Fermi level on the right evolves to be at R at r again, but its left-behind complement 
(a negative charge) has evolved to be at L. (Tlr{t) in this case means the transition of an 
empty state or positively charged particle from right to left (Fig. 5b). Although one docs not 
usually introduce the concept of holes with transport in normal conducting metals, it makes 
sense here to call unoccupied states simply "holes". In this way the model already includes 
the dual nature of charge carriers needed for the superconducting case. In the normal 
conducting case states do not change in nature (electron or hole) or energy during their 
evolution (little zigzags are drawn in Fig. 5 only to make the multiple hoppings visible). At 



I/e 



= (^c,LR (1 + G'^j^L^lr) 9rr O-rlGI^ 

+ cfc,lr G^'jij^aRL gtl (1 + ctlrG^l) 

+ cTcRL Gl^aLR Qrr (1 + ctrlGIr) 

+ (^c,RL (1 + Glj^aLR) qIl CTlrGrr 



(43) 
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the superconducting junction, Andreev reflection can be interpreted as changing an electron 
into a hole or vice versa and mirroring its energy at the Fermi level ^\ (see Fig. 9 in section 
6). To have a charge carrier at a certain energy level at time r to make a certain transition, 
it is important that there was one at the corresponding energy in the original uncoupled 
system. The coupling may have changed the distribution compared to the occupation in 
uncoupled bulk reservoirs. And the applied voltage imposes a non-equilibrium situation, 
anyway. For the evolution of a state from the right as shown in Fig. 5 it does not matter 
whether the state at the respective energy on the left is occupied or empty. If, for example 
the regarded energy lies below the Fermi levels both left and right, there will be two states 
evolving as an electron on the plus and a hole on the minus branch of the Keldysh contor, 
one having originated at R and the other at L. These original, uncoupled and independent 
states are our basis, especially for calculating an expectation value as trace. They do not 
interfere. Terms with g^J^ and g~l^£ are simply added in (43). Schemes analogous to Fig.5 
could be drawn for the terms from the last three lines of (43) as well. The conclusion of the 
whole argumentation of how to let the original Fermi occupation function for the reservoirs 
left and right enter the current calculation is that g~^£ has to change sign at the left Fermi 
energy and g^j^ at the right Fermi energy. Let us note g^~ like g^' and g"" for any bulk 
reservoir with Fermi level at u = 0. We shall keep g^^ corresponding to g^ and g"' as in 
(40) for occupied electron states below the Fermi level and change the sign for empty states 
above it. 

J) r ■ , „ i for Co" < 

g+-{r, r') = j du g^- {u) e--^^^^ ) with g+-{u) = I (44) 

— i for c<j > 

Although it might be practical to use (19), (23), (44), (18) and (36) in (41) to quite directly 
produce an expression that calculates the current finally as an integral over frequency and 
in our simple model can even be analytically evaluated, in parallel to |58] we shall use the 
transfer Green's functions here. The current is best translated into an expression of the 
transfer functions from the form already resolved into LR-components (43). Furthermore 
eliminate o"c through o"c,l_r = —ctlr and a^Ri = ctrl- L- and R-indices follow logically from 
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(26), for example Gllctlr = QllTlr- 

I/e = a^R (1 + g'nnTl,^) g^^ T^lQIl 

+ (^LRgWT'RLatli'^ + TlM 

- ^Rl91lTIh9rr ii + T^RR9RR) 

- URL (1 + gl^Tl^) gtE TUg^RR (43a) 

As quantities here are no longer matrices, but simply functions of time or frequency, the 
leading and in the second and third line of (43a) can be moved to the ends of 
the products. Prom (41) we remember that their argument is r, the same as the second 
argument of the last factor to the right. Then using relations from (27a) and complementary 
forms the current formula simplifies to 

^/^ — '^LRgRR'^RigLL + gRR^RLgLL '^LR 

- 9lLTLRgtRTRL - T'kL9tLTlR9RR (43b) 

The terms with the factors 1 from the brackets have elegantly been made to vanish. Now 
we use the Fourier representations for all functions in (43b). As all terms follow the same 
scheme, the second is treated in an exemplary way here (all integrals run from —oo to +cxd): 



j dri dT2 dT3 dwi duj duj2 duj3 j^gRR{uJi)e "^^'^ e^'^'^'^'^ -^:^T^l{'^2) 

Doing the Ti,T2- and T3-integrals gives {2'k)^5{u!i —lo^ — eV/h)S{u!2 — 00)6(00 — cus). Then even 
the exponentials with r cancel and the term simplifies to 

1 

2^ 



du; g^^^iu; + eV/h) 9tM TU^) 

In our case (31) tells us that T^R/RLii^) are all identical and real, g^L/RRioo) — —i is the 
complex conjugate of g^L/RRi^) ~ ^ ~9ll/rr complex conjugate of g^L/RR the 

same lo as g^^ {uo) is purely imaginary, too. Thus the last two terms in (43b) are the complex 
conjugates of the first two and thus twice the real part of these first two can be taken for 
I/e. In the superconducting version of the model, where T and g actually are a;-dependent. 
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complex conjugate relations [59] also exist between Ts as well as gs, and the current formula 
can be reduced in the same way. Just to note the quite general formula |^ in short form: 

I/e = 2 Re {Tl^g^-T^Lgl^ + g'nBT'nLgtlTU} (45) 

For the integrand from the second term in our model we get 

g^j,^{u + eV/h)TJ,,ico)gtE{co)Tlniu^) = = T- 



1+^2^ ' 'l+t2 ' (1+^2)2 

where the signs refer to uj greater or less than zero and come out reversed for the first 
term with g^]^, because there is g"" instead of g^ . Care has to be taken with the reference 
point for uj in both terms. This may easily be overlooked in the normal conducting case 
here in contrast to the superconducting case where T indeed is a;-dependent and like G as a 
function of only one argument always referred to the same Fermi level (the left, for example). 
If we call the argument of (7^^ from the second term uj, the one for g^^ in the first term is 
UJ + eV/h. On an C(j(i)-axis, the second term changes sign at zero, however, the first jumps 
at —eV/h (Fig.6). A shift of the integration parameter cannot be made independently for 
both terms. Thus, for the normal conducting model here 

I/e = 2Re^j duj{ TUuj)g^-iuj + ^^)nd^)gl^i^) 

+ + ''^)T'aL{^)9tMTUi^) } (45a) 

In principle the convention is needed, that the T-argument always refers to the left side, 
but for g general formula like (44) for a single bulk with Fermi level at zero frequency are 
applied. A more involved situation where choosing integration intervals consistently for all 
contributing current terms is crucial can be found in [59] . 

From Fig.6 it is easily seen that that the integral is twice the constant f^i^f^yi integrated 
over an interval of length eV/h, and outside that interval contributions cancel. With the 
other factor 2 from twice the real part and the prefactor l/27r the result for the current 
finally is 

^/"^ " ^ ■ 2^ ■ ^ ■ (1 + t2)2 ■ ~n ^ 

V_h_ ( 4t2 \ ]^_ e2 4t2 

The factor 9 = 4t2/(l + 12~)2 ^j^g conductance in units of the quantum conductance or its 



inverse the resistance in units of the quantum resistance [58] • (The conductance of a channel 
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FIG. 8: (a) Evolution of an electron state from R such that at r the charge is transferred to the 
left, (b) Evolution of an unoccupied state above the Fermi level from R such that at r it gets filled 
by a charge from the left or the hole state shifts to the left. The lower pictures sketch the effect 
of the original uncoupled state for the current. Sending electrons from the right to the left will be 
compensated by electrons sent from left to right. But with levels as in Fig.lc for holes sent from 
right to left there will be a range between the two Fermi levels where there are no counterbalancing 
holes going the other way. (For the left, electrons going to the right are not outweighed in this 
range.) 

51 

















eV/Pr 


(0 


^ '■ ► 



FIG. 9: The two terms under the integral from equation (45). They cancel except on an interval 
of length eV/ h. 

doubles if two spin states are allowed. Then h/2e^ — 13/cQ should be taken as resistance 
unit.) 6 is the transmission probability of the conductance channel through the junction 
we regarded. The result that in the normal-conducting case the current is proportional to 
the voltage is not at all surprising, of course. The non-trivial result is the conversion of the 
quantum mechanical transmission amplitude t to the measurable transmission probability 
9. 9 — Q lit = And a totally open channel with t — 1 has transmission probability 9 — 1. 

It may seem a contradiction on the one hand calculating the current from a changing 
amount of charge on one side and on the other hand saying that missing charges are replaced 
and superfluous ones led away by the voltage source. A slightly different viewpoint may help 
to get convinced that the calculated quantity is indeed the current in the stationary, but 
non-equilibrium system. The crucial point was putting on to evaluate dp^/dt in state ^E'(t). 
Without further ado we could not tell whether this state of the coupled system was occupied 
or not. The Schrodinger equation (5) set up the left and right material properties as well 
as the coupling across the junction, however, did not take any account of the effect of the 
voltage source. Without need to specify real locations for the division, just principally view 
our structure as consisting of a junction region and leads. The ^(t) defined through (5) 
describes states in the junction region. But think of them as offered by the system and 
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following their time development whether occupied or not. The leads are always occupied 
exactly up to their Fermi levels. Electrons freshly supplied by the voltage source need 
not be in phase with present ones. A random phase is most easily modelled by assuming 

the left lead wave function at any time without a phase in contrast to the assumed 

exp{—iujLT) \ 

for the left side of the junction region. Which states to which extent actually 

° J 

get occupied in the left side of the junction region, that is piit), is determined by the overlap 
of \t'(t) with the phaseless left lead wave function ( | , which we recognize as \t'5,(— C)o). 

W 

Pl = \ < i III ^''M>p. The current then is the change in time of pi due to charge 

flow through the junction, that is processes inside the junction region only, corresponding 
to what our Hamiltonian was set up for. A similar line of thought with state overlaps can 
be applied to the passing of charges out of the junction region into the (right) lead; it may 
be helpful to view this as putting empty states or holes into the junction region, though. Of 

course, only the L-part of \l'(t) can overlap with . can thus be replaced by 

'^{t) in our new expression for p^. Recalling the operator definition (34) the time derivative 
of pL 

i--ii<(;)i*w>i^-i{<*<*)i(;)>< (;)i*w> 

is identical to the ansatz made by putting together (33) and (34) in (37). 



III. 6. Comment on the two-level system 

At first glance the problem posed by the Schrodinger equations without and with coupling, 
(1) and (5), especially if we regard a single energy level on each side with corresponding 
ujl and ujr as in Fig.lc, looks like the two-level system known from quantum mechanics 
textbooks [HI]- By the coupling the two energy levels are shifted. The new eigenstates lie 
further apart from one another. If the system is initiated in one of the uncoupled states, it 
will oscillate harmonically between the two levels, and both both the period of the oscillation 
as well as the maximum transition probability to the other state depend on the energy 
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FIG. 10: Hypothetical division of the structure into the actual junction region and bulk leads. 

difference of the original states. 

Even if the time dependence of the interaction (18) can be got rid of by changing to the 
interaction picture, it makes no sense to numbly calculate eigenvalues and eigenvectors of 

a Hamiltonian . It is not clear from which reference levels such eigenenergies 

E± should be counted as zero levels are different for lol and lor. If one tries to refer the left 



eV 
2 



and right energies already of the uncoupled system to the same reference level, Ef{R) + 
for example (Fig.8), (and ^5?) cannot be taken as eigenstates any more, because from 
ihdre-^'^i^^ = hwLe-''^^^ with cul = cj^ + ^ it does not follow that ihdre-''^<+^^/^^^ equals 
^/^g-j(aj^+ey/2)T r^Yie standard treatment of the two-level system is not applicable to the 
non-equilibrium situation. 



III. 7. Superconducting junction 

Despite having been the simplest example to introduce the Green's functions scheme, 
applying the formalism to the normal conducting junction to get out that the current is 
proportional to the applied voltage was breaking a butterfly on the wheel. Although repeat- 
edly mentioned, fully developing the extension to the superconducting case is beyond the 
scope of this presentation. But the resulting current-voltage characteristics shall be shown 
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FIG. 11: Attempt to adjust energy levels from left and right to same reference. 



as a plea for the usefulness of the method. A different approach based on matching wave 
functions [63] leads to identical results, though. 

Some formula shall be listed, because they are not necessarily written out in complete 
analogy to the presentation here in [SS] and other literature. Working in the quasiparticle 
picture, in the superconducting case each entry in LR-space of a Green's or transfer function 
expands into another 2x2 matrix in Nambu space over electrons and holes. There is 



rla I \ 



gee geh 
ghe ghh 



r/a 



[UJ 



-uj^irj A/h 
—A/h uj±iri 



(47) 



LL/RR 

g is the same for LL and RR. Different signs refer to the retarded and advanced function. 
A is half the gap of the superconductor. The small imaginary part rj is added to get the 
correct root besides slightly smoothening singularities. g{T,T') = V/h J duj g{uj) e"*"^"^ e*'^'^' 



still holds. 



C" LR/RL 



a 



hh 



LR/RL 



h 1 1 e^'^^^/^ 
V 1 _^^TieVr/h 



(48) 
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FIG. 12: Illustration of Andreev reflection as mirroring at the Fermi levels for electron-hole con- 
version and vice versa. For either side the incoming and reflected levels can lie inside or outside 
the gap. Only the beginning and end of a complete multiple-reflection process have to be in the 
electron or hole reservoir below or above the gap, respectively. 

Here signs refer to the LR- and RL-direction, respectively. Remark the reversed signs for 
holes in the exponential compared to electrons. There is no electron-hole conversion during 
a single hopping. The eh- and /le-elements of g introduce so-called Andreev reflection. 
(Cooper-pair tunneling [65] is not included.) T is a full matrix in e/i-space and has to be 
set up as 

T) 1 „ /• _ 

(49) 



ineVr/h ^^t' 



We skipped the distinction of reference ujjik as in (25), because in practice only one type, 
e.g. Ti/j, has to be calculated, as then Tr/, follows from complex conjugation, n runs over 
all integers. A Dyson equation like (32) leads to a recursion which connects T„ to T„_2 and 
T„+2- It can be solved by reasonably truncating the n-range. Then, just also summing over 
n, the established current formula (45) can be used. 

Expressions like (32) have to sum over all possible combinations of e- and /i-indices. For 
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FIG. 13: Calculated current-voltage curves for transport channels of three different transmissions 
(normalized to 9) [reproduced after [58]]. The normal conducting IV is added for comparison (keep 
A just for units). 



example: 



rjiee ee i ee ee ee ee /jiee , ee eh hh he ijiee 

-'-LR — '-'LR "T LRiJ RR'^ RLiJ LL-'- LR "r LrS RR'^ RlS LL-'- LR 



I „ee ^ee „ee „eh rphe i „ee „eh „hh ^hhrphe (Kr\\ 

+ (^lr9rr(^rl9lL'^lr + (^lr9rr(^rl9lL'^lr (50) 

To point out the Andreev reflection in the model, we look at a term gagcg with alternating 
e- and /i-indices from the development (21). Such a sequence will also be contained as parts 
in higher-order terms from (21) or (24). Contracting inner time arguments results in 

dr, dT2 du2 duj^ du jgfR{uj2)e~'^'^e'^^^' _t^-ieVr,/h . 
'^gf^Me-'-^^^e'-^^- l(_t)e-^-/^ |^M^)(^)e— ^e^-' 

= (—t ) — I duj duoi duj2 o{uj2 uJi) o{uJi u) ■ 

h J h n 

9^Ri'M gfLM g'i^^uj) e-' 

= -t'jl du gfn{u + 2^) gtd^ + f) gf^{u) e-^(^+-^/^)^ e^' 

(51) 

because the ri and r2-integration give 47r^(5(co'2 — eV/h — uji)5{ijJi — eV/h — oj). In Fig.9 for 
holes the energy axis has to be reversed as compared to electrons. Therefore the initial uj is 
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negative. For the argument here, it does not matter if the beginning factor in (51) is g'^j^ or 
g^j^. If it were h, then the Andreev reflection marked by the dashed arrow would follow. The 
last index of the last g^R could also be h instead of e. Whereas in the normal conducting 
case in multiple reflections the charge carrying particle always came back to the same energy 
with regard to the Fermi level when coming back to the same side of the junction (Fig.l), in 
(51) the frequency argument of the first gnR is shifted by twice the applied voltage equivalent 
with respect to the second g^ji (Fig.9). Such shifts enable transport even when the gaps of 
the left and right side are still overlapping, that is for voltages below 2A/e (Fig. 10). In the 
Jl^-curves steps towards lower voltages at fractions of 2A are associated to ever higher-order 
Andreev reflections. Of course, these current contributions become the more prominent the 
greater the channel transmission. Besides the step heights and positions the curvature on 
each step is an important signature. The model agrees excellently with experimental results 
[53] . The interference of transport processes of different order is correctly taken into account. 
Simpler models fB6j that add up tranfer probabilities proportional to 6°'"'^'^'^ cannot reproduce 
these characteristics and are only valid for low transmissions << 1. 



IV. CONCLUSIONS AND OUTLOOK FOR SECTION III 

It has been explained on a quite basic level how quantum transport between two reservoirs 
in a stationary non-equilibrium state can be modelled. The purpose was to present a Green's 
functions technique for handling coupling in the context of a field of current research interest, 
namely transport through point contacts. Although besides the general formalism (also 
see section I) requiring the development of quite some more subject specific mathematical 
framework such as the transfer Green's functions, the transport through a contact with 
arbitrary transmission is a suitable example to illustrate the inclusion of interaction up to 
all orders in the implicit Dyson equations. Calculations have been carried through in every 
detail for the normal-conducting single junction. Basic formula and results have been given 
for the superconducting junction. The decomposition into transport channels (eigenmodes) 
of a point contact can be inferred from its superconducting transport characteristics which 
can be taken like a PIN-code [M]- The presented Green's functions formalism has a great 
potential for extension. Systems of two more or less coherently linked junctions can be 
modelled [59]. Models for transport through molecules [671 l68l [69] or atomic chains [70] so 
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far mainly rely on ab initio calculations of the density of states. Time dependent density 
functional theory for non-equilibrium situations is also developed j71j . 
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